Children - Study 2 - Results

This is the results section for the Study 2 NSE & SE CHILDREN watching ASL Stories. We have two main factors:

1. Language (Sign v. English)
2. Direction (Forward v. Reversed)

We are taking out one older KODA (Ethan, 10.5 yrs) to balance the groups better.

Demographics

library(tidyverse)
library(janitor)
library(lme4)
library(lmerTest)
library(scales)
library(feather)
library(GGally)

kids <- read_feather("cleanedchildeyedata.feather") %>%
#  mutate(age = age*12) %>%
  select(participant, language, age, gender, story, direction, mark, trial, repetition, aoi, secs, percent) %>%
  rename(name = participant) %>%
  filter(age < 9) %>% # Take out Ethan
  # mutate(agegroup = case_when(
  #   age <= 8.99 ~ "younger",
  #   age >= 9.0 & age < 15 ~ "older"
  # )) %>%
  # filter(!is.na(agegroup)) %>%
  mutate(language = case_when(
    language == "english" ~ "NSE",
    language =="sign" ~ "SE"
  )) %>%
  rename(lang = language)

kidsinfo <- kids %>%
  select(name, lang, age, gender) %>%
  distinct() %>%
  group_by(lang) %>%
  summarise(N = n(),
            age_mean = mean(age),
            sd = sd(age),
            min = min(age),
            max = max(age))

genders <- kids %>%
  select(name, lang, age, gender) %>%
  distinct() %>%
  group_by(lang, gender) %>%
  summarise(N = n()) %>%
  spread(gender, N)

kidsinfo <- left_join(kidsinfo, genders) %>%
  select(lang, N, Female, Male, age_mean, sd, min, max) %>%
  print()

# babies$agegroup <- fct_relevel(babies$agegroup, c("younger","older"))


# IF we do age groups, use this code
# 
# babiesinfo <- babies %>%
#   select(name, lang, age, agegroup, gender) %>%
#   distinct() %>%
#   group_by(lang, agegroup) %>%
#   summarise(N = n(),
#             age_mean = mean(age),
#             sd = sd(age),
#             min = min(age),
#             max = max(age))
# 
# genders <- babies %>%
#   select(name, lang, age, agegroup, gender) %>%
#   distinct() %>%
#   group_by(lang, agegroup, gender) %>%
#   summarise(N = n()) %>%
#   spread(gender, N)
# 
# babiesinfo <- left_join(babiesinfo, genders) %>%
#   select(lang, agegroup, N, Female, Male, age_mean, sd, min, max) %>%
#   print()

Let’s plot the ages, and check if there is significant difference in ages between the two groups?

# Boxplot
kids %>%
  select(name, age, lang) %>%
  distinct() %>%
  ggplot(aes(x = lang, y = age, fill = lang)) + geom_boxplot(width = 0.5) + guides(fill = FALSE)

kids %>%
  select(name, age, lang) %>%
  distinct() %>%
  ggplot(aes(x = age, fill = lang)) + geom_histogram() + facet_grid(lang ~ .)

# T-test
nse_age <- kids %>% filter(lang == "NSE") %>% select(name, age) %>% distinct()
se_age <- kids %>% filter(lang == "SE") %>% select(name, age) %>% distinct()
t.test(nse_age$age, se_age$age)

    Welch Two Sample t-test

data:  nse_age$age and se_age$age
t = 0.14316, df = 32.568, p-value = 0.887
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.020020  1.174354
sample estimates:
mean of x mean of y 
 4.986667  4.909500 

Global Looking

For children, we calculated percentages based on overall clip length as the denominator. In this way, we can meaningfully contrast looking times at the videos (which are variable lengths) based on different factors. But when we go to AOI analysis we need to re-calculate the percentages so the denominator is based on total looking time, not overall clip length.

The chart below shows there seems to be an effect of age; older kids look longer at it than younger kids. Maybe not too surprising. It means we need to keep age in any models we run. Let’s analyze a bit more below.

kids$lang <- as.factor(kids$lang)
kids_overall_looking <- kids %>%
  group_by(name, age, lang, direction, story, repetition) %>%
  summarise(percent = sum(percent)) # gets total looking percent for each trial for each kid

# Table of means
kids_overall_looking %>% 
  group_by(name, lang, direction) %>%
  summarise(percent = mean(percent)) %>% # get average looking percent for each kid
  group_by(lang, direction) %>%
  summarise(mean_percent = mean(percent),
            count = n(),
            sd = sd(percent),
            se = sd/sqrt(count)) %>%
  select(-sd) %>%
  print()

ggplot(kids_overall_looking, aes(x = age, y = percent, color = direction, fill = direction)) + 
  geom_jitter(alpha = 0.5) +
  facet_grid(. ~ lang) +
  geom_smooth(method = "lm", se = TRUE) +
  ggtitle("Video Attention") +
  xlab("age (months)") +
  ylab("percent looking") + 
  theme_bw() + 
  scale_y_continuous(limits = c(0,1), labels = percent) 

# Plot
# babies_overall_looking %>% 
#   group_by(lang, direction, name) %>%
#   summarise(percent = mean(percent)) %>% # gets average looking percent for each baby
#   group_by(lang, direction) %>%
#   summarise(mean_percent = mean(percent), # gets group averages
#             count = n(),
#             sd = sd(percent),
#             se = sd/sqrt(count)) %>% 
#   ggplot(aes(x = lang, y = mean_percent, fill = direction)) + 
#   geom_col(position = "dodge") + 
#   geom_errorbar(aes(ymin = mean_percent - se, ymax = mean_percent + se), 
#                 position = position_dodge(width = 0.9), width = 0.25) + 
#   scale_y_continuous(limits = c(0,1), labels = percent) +
#   theme_minimal() + 
#   theme(panel.grid.major.x = element_blank()) +
# #  facet_wrap("lang") +
#   ggtitle("Video Attention") +
#   xlab("") +
#   ylab("percent looking")

# babies_overall_looking %>%
#   ggplot(aes(x = lang, y = percent, fill = direction)) +
#   facet_wrap("agegroup") + 
#   geom_violin()

A linear model shows a significant effect of age. Overall, Age seems to increase overall looking by about 3% every year. However, there are no differences between NSE v. SE, or reversal, on how long they looked, so that’s good.

global_lm <- lmer(percent ~ age + lang * direction + (1|name) + (1|story), data = kids_overall_looking)
summary(global_lm)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: percent ~ age + lang * direction + (1 | name) + (1 | story)
   Data: kids_overall_looking

REML criterion at convergence: -112.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.4984 -0.6942  0.2184  0.7802  2.5041 

Random effects:
 Groups   Name        Variance  Std.Dev.
 name     (Intercept) 0.0071122 0.08433 
 story    (Intercept) 0.0003727 0.01931 
 Residual             0.0397608 0.19940 
Number of obs: 471, groups:  name, 35; story, 8

Fixed effects:
                           Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)                0.582019   0.062007  34.950676   9.386 4.39e-11 ***
age                        0.032069   0.010837  30.892721   2.959  0.00588 ** 
langSE                     0.032272   0.039174  51.060608   0.824  0.41388    
directionreversed         -0.031148   0.028318 375.551029  -1.100  0.27206    
langSE:directionreversed  -0.003914   0.037809 384.555021  -0.104  0.91760    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) age    langSE drctnr
age         -0.872                     
langSE      -0.376  0.019              
dirctnrvrsd -0.233  0.004  0.368       
lngSE:drctn  0.179 -0.006 -0.481 -0.756

confint(global_lm)

Computing profile confidence intervals ...

                               2.5 %     97.5 %
.sig01                    0.05357505 0.11195261
.sig02                    0.00000000 0.05071616
.sigma                    0.18641094 0.21320556
(Intercept)               0.46258809 0.70154684
age                       0.01118730 0.05298192
langSE                   -0.04342330 0.10755115
directionreversed        -0.08722835 0.02442038
langSE:directionreversed -0.07771906 0.07088373

#ggcoef(global_lm)

library(sjPlot)

Learn more about sjPlot with 'browseVignettes("sjPlot")'.

library(sjmisc)

package ‘sjmisc’ was built under R version 3.6.2Learn more about sjmisc with 'browseVignettes("sjmisc")'.

Attaching package: ‘sjmisc’

The following objects are masked from ‘package:janitor’:

    remove_empty_cols, remove_empty_rows

The following object is masked from ‘package:purrr’:

    is_empty

The following object is masked from ‘package:tidyr’:

    replace_na

The following object is masked from ‘package:tibble’:

    add_case

library(sjlabelled)

Attaching package: ‘sjlabelled’

The following object is masked from ‘package:forcats’:

    as_factor

The following object is masked from ‘package:dplyr’:

    as_label

tab_model(global_lm, show.se = T, show.stat = T)

AOI Looking

Now we’ll re-calculate the percentages so the denominator is based on total looking time. All AOIs should sum up to 100% for each trial and each baby. Next let’s make a boxplot of all AOIs. Interesting, definitely more MidFaceBottom focus here than we had with babies, but also more distribution too.

# Recalculate percent
kids <- kids %>% 
  ungroup() %>%
  select(-percent) %>%
  group_by(name, lang, age, direction, story, mark, trial, repetition, gender) %>%
  mutate(totalsec = sum(secs)) %>%
  group_by(name, lang, age, direction, story, mark, trial, repetition, gender, aoi) %>%
  summarise(percent = secs/totalsec)

# Boxplot
kids %>%
  ggplot(aes(x = aoi, y = percent, fill = direction)) + 
  geom_boxplot() +
  ggtitle("AOI Attention") +
  theme_bw() + 
  xlab("") +
  theme(axis.text.x = element_text(angle=45, hjust = 1),
        panel.grid.major.x = element_blank()) +
  scale_y_continuous(labels = scales::percent, limits = c(0,1))

It appears two important AOIs are MidChestTop and MidFaceBottom. Let’s look again only at midline AOIs:

midline = c("Belly","BelowChest","MidChestBottom","MidChestCenter","MidChestTop",
            "MidFaceBottom","MidFaceCenter","MidFaceTop")
kids %>%
  filter(aoi %in% midline) %>%
  ggplot(aes(x = aoi, y = percent, fill = direction)) + 
  geom_boxplot() +
  ggtitle("Midline AOI Attention") +
  theme_bw() + 
  xlab("") +
  theme(axis.text.x = element_text(angle=45, hjust = 1),
        panel.grid.major.x = element_blank()) +
  scale_y_continuous(labels = scales::percent, limits = c(0,1))

I’m going to run linear models with only MidChestTop or MidFaceBottom, and see what happens. No age interactions.

MidChestTop:

• No effect of age.
• No effect of language.
• Weak effect of direction (p = 0.084) - reversed means they look ~ 4.7% less at midchesttop.
• No language X direction interaction.

MidFaceBottom:

• No effect of age.
• Significant effect of language (p = 0.004) - SE look at MidFaceBottom +19% more than NSE children
• No effect of direction.
• No language X direction interaction.

kids %>%
  filter(aoi %in% c("MidFaceBottom","MidChestTop")) %>%
  ggplot(aes(x = age, y = percent, color = direction, fill = direction)) + 
  geom_jitter(alpha = 0.5) + 
  geom_smooth(method = "lm", se = FALSE) +
  scale_y_continuous(limits = c(0,1), labels = percent) +
  theme_bw() + 
#  theme(panel.grid.major.x = element_blank()) +
  facet_grid(aoi ~ lang) +
  ggtitle("AOI Attention") +
  xlab("") +
  ylab("percent looking")

midchesttop_lm <- lmer(percent ~ age + lang * direction + (1|name) + (1|story), data = filter(kids, aoi == "MidChestTop"))
summary(midchesttop_lm)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: percent ~ age + lang * direction + (1 | name) + (1 | story)
   Data: filter(kids, aoi == "MidChestTop")

REML criterion at convergence: -108

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.0759 -0.6139 -0.1635  0.5257  3.4431 

Random effects:
 Groups   Name        Variance  Std.Dev.
 name     (Intercept) 0.0236542 0.15380 
 story    (Intercept) 0.0001076 0.01037 
 Residual             0.0375728 0.19384 
Number of obs: 471, groups:  name, 35; story, 8

Fixed effects:
                           Estimate Std. Error         df t value Pr(>|t|)   
(Intercept)                0.332735   0.097932  33.690232   3.398  0.00176 **
age                       -0.004869   0.017508  32.325032  -0.278  0.78269   
langSE                    -0.067047   0.058495  39.320582  -1.146  0.25864   
directionreversed         -0.047327   0.027264 337.976644  -1.736  0.08350 . 
langSE:directionreversed   0.049457   0.036436 348.971721   1.357  0.17554   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) age    langSE drctnr
age         -0.892                     
langSE      -0.360  0.022              
dirctnrvrsd -0.142  0.003  0.235       
lngSE:drctn  0.109 -0.004 -0.310 -0.751

confint(midchesttop_lm)

Computing profile confidence intervals ...

                               2.5 %      97.5 %
.sig01                    0.11354190 0.193465409
.sig02                    0.00000000 0.042445622
.sigma                    0.18117827 0.207111281
(Intercept)               0.14404353 0.521394766
age                      -0.03861109 0.028872872
langSE                   -0.17970025 0.045555254
directionreversed        -0.10070179 0.006183759
langSE:directionreversed -0.02216830 0.120665650

#ggcoef(midchesttop_lm)

midfacebottom_lm <- lmer(percent ~ age + lang * direction + (1|name) + (1|story), data = filter(kids, aoi == "MidFaceBottom"))
summary(midfacebottom_lm)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: percent ~ age + lang * direction + (1 | name) + (1 | story)
   Data: filter(kids, aoi == "MidFaceBottom")

REML criterion at convergence: 27

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.69174 -0.65300 -0.01245  0.70728  2.54595 

Random effects:
 Groups   Name        Variance Std.Dev.
 name     (Intercept) 0.025609 0.16003 
 story    (Intercept) 0.001335 0.03654 
 Residual             0.050261 0.22419 
Number of obs: 471, groups:  name, 35; story, 8

Fixed effects:
                           Estimate Std. Error         df t value Pr(>|t|)   
(Intercept)                0.359105   0.104204  34.706091   3.446  0.00151 **
age                       -0.005517   0.018452  32.182394  -0.299  0.76685   
langSE                     0.188370   0.062416  41.098544   3.018  0.00436 **
directionreversed         -0.005449   0.032290 418.178679  -0.169  0.86607   
langSE:directionreversed  -0.049419   0.043090 421.686350  -1.147  0.25208   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) age    langSE drctnr
age         -0.883                     
langSE      -0.361  0.021              
dirctnrvrsd -0.158  0.002  0.266       
lngSE:drctn  0.122 -0.003 -0.345 -0.762

confint(midfacebottom_lm)

Computing profile confidence intervals ...

                               2.5 %     97.5 %
.sig01                    0.11698914 0.20280761
.sig02                    0.00000000 0.07849412
.sigma                    0.20953296 0.23957130
(Intercept)               0.15817434 0.56034402
age                      -0.04117347 0.03009216
langSE                    0.06791089 0.30869148
directionreversed        -0.06926187 0.05753800
langSE:directionreversed -0.13345107 0.03574359

#ggcoef(midfacebottom_lm)

# Bar chart
# babies %>%
#   filter(aoi %in% c("MidFaceBottom","MidChestTop")) %>%
#   group_by(agegroup, lang, direction, name, aoi) %>%
#   summarise(percent = mean(percent)) %>% # gets average looking percent for each baby
#   group_by(agegroup, lang, direction, aoi) %>%
#   summarise(mean_percent = mean(percent), # gets group averages
#             count = n(),
#             sd = sd(percent),
#             se = sd/sqrt(count)) %>% 
#   ggplot(aes(x = lang, y = mean_percent, fill = direction)) + 
#   geom_col(position = "dodge") + 
#   geom_errorbar(aes(ymin = mean_percent - se, ymax = mean_percent + se), 
#                 position = position_dodge(width = 0.9), width = 0.25) + 
#   scale_y_continuous(limits = c(0,1), labels = percent) +
#   theme_minimal() + 
#   theme(panel.grid.major.x = element_blank()) +
#   facet_grid(aoi ~ agegroup) +
#   ggtitle("Video Attention") +
#   xlab("") +
#   ylab("percent looking")

Face-Chest Ratio

Next, we’ll define a Face-Chest Ratio (FCR) such that:

1. MidFaceCenter, MidFaceBottom = Face
2. MidChestTop, MidChestCenter, MidChestBottom, BelowChest = Chest
3. FCR = face - chest / face + chest

We did not include Belly or MidFaceTop because of very low looking rates according to the boxplots above.

kids_fcr <- kids %>%
  ungroup() %>%
  spread(aoi,percent) %>%
  group_by(name, age, lang, gender, direction, story, repetition) %>%
  summarise(face = sum(MidFaceCenter, MidFaceBottom, na.rm = TRUE),
         chest = sum(MidChestTop, MidChestCenter, MidChestBottom, BelowChest, na.rm = TRUE),
         fcr = (face - chest) / (face + chest))

# Table of means
kids_fcr %>% 
  group_by(lang, direction, name) %>%
  summarise(fcr = mean(fcr)) %>% # gets average looking percent for each baby
  group_by(lang, direction) %>%
  summarise(mean_fcr = mean(fcr), # gets group averages
            count = n(),
            sd = sd(fcr),
            se = sd/sqrt(count)) %>%
  select(-sd) %>%
  print()

kids_fcr %>% 
  group_by(lang, name) %>%
  summarise(fcr = mean(fcr)) %>% # gets average looking percent for each baby
  group_by(lang) %>%
  summarise(mean_fcr = mean(fcr), # gets group averages
            count = n(),
            sd = sd(fcr),
            se = sd/sqrt(count)) %>%
  select(-sd) %>%
  print()

# Plot
ggplot(kids_fcr, aes(x = age, y = fcr, color = direction, fill = direction)) + 
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_jitter(alpha = 0.5) + 
  geom_smooth(method = "lm", se = FALSE) +
  scale_y_continuous(limits = c(-1,1)) +
  theme_bw() + 
  theme(panel.grid.major.x = element_blank()) +
  facet_grid(. ~ lang) +
  ggtitle("Face-Chest Ratios") +
  xlab("") +
  ylab("FCR")

# Bar chart
# babies_fcr %>% 
#   group_by(agegroup, lang, direction, name) %>%
#   summarise(fcr = mean(fcr)) %>% # gets average looking percent for each baby
#   group_by(agegroup, lang, direction) %>%
#   summarise(mean_fcr = mean(fcr), # gets group averages
#             count = n(),
#             sd = sd(fcr),
#             se = sd/sqrt(count)) %>% 
#   ggplot(aes(x = lang, y = mean_fcr, fill = direction)) + 
#   geom_col(position = "dodge") + 
#   geom_errorbar(aes(ymin = mean_fcr - se, ymax = mean_fcr + se), 
#                 position = position_dodge(width = 0.9), width = 0.25) + 
#   scale_y_continuous(limits = c(-1,1)) +
#   theme_minimal() + 
#   theme(panel.grid.major.x = element_blank()) +
#   facet_wrap("agegroup") +
#   ggtitle("Face-Chest Ratios") +
#   xlab("") +
#   ylab("FCR")

What will a linear mixed model tell us? (with no age interactions)

• No effect of age. Interesting. Maybe just becuase we don’t have that many babies.
• Effect of language: SE babies have overall +0.36 FCR than NSE babies (p = 0.03). SE looks at the face more than the chest.
• No effect of direction. Interesting.
• No interaction. Interesting.

fcr_lm <- lmer(fcr ~ age + lang * direction + (1|name) + (1|story), data = kids_fcr)
summary(fcr_lm)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: fcr ~ age + lang * direction + (1 | name) + (1 | story)
   Data: kids_fcr

REML criterion at convergence: 747.9

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-3.04087 -0.55798  0.04156  0.68150  2.93235 

Random effects:
 Groups   Name        Variance Std.Dev.
 name     (Intercept) 0.17809  0.4220  
 story    (Intercept) 0.01392  0.1180  
 Residual             0.22797  0.4775  
Number of obs: 471, groups:  name, 35; story, 8

Fixed effects:
                           Estimate Std. Error         df t value Pr(>|t|)  
(Intercept)               -0.051211   0.268645  34.803206  -0.191   0.8499  
age                       -0.001861   0.047533  32.187297  -0.039   0.9690  
langSE                     0.357443   0.157942  38.317004   2.263   0.0294 *
directionreversed         -0.047553   0.069424 433.299130  -0.685   0.4937  
langSE:directionreversed  -0.068533   0.092601 434.179164  -0.740   0.4596  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) age    langSE drctnr
age         -0.882                     
langSE      -0.355  0.022              
dirctnrvrsd -0.132  0.001  0.227       
lngSE:drctn  0.102 -0.002 -0.293 -0.766

confint(fcr_lm)

Computing profile confidence intervals ...

                               2.5 %     97.5 %
.sig01                    0.31504728 0.53086491
.sig02                    0.05652398 0.22952050
.sigma                    0.44621934 0.51015034
(Intercept)              -0.56992177 0.46783485
age                      -0.09379276 0.09001524
langSE                    0.05243219 0.66250287
directionreversed        -0.18416110 0.08799143
langSE:directionreversed -0.24911301 0.11393372

#ggcoef(fcr_lm)

write_csv(kids_fcr, "fcr_trial_level_values_children.csv")

tab_model(fcr_lm, show.se = T, show.stat = T)

The Large Table

I would like a large table with all individual percent looking means for each AOI and the individual FCR values, with ages, gender, video group for each child. (collapsed across stories and trials)

# Collapse across stories and trials 

kids_spread <- kids %>%
  group_by(name, lang, age, gender, direction, aoi) %>%
  summarise(percent = mean(percent, na.rm = T)) %>%
  spread(aoi, percent)

kids_fcr_spread <- kids_fcr %>%
  group_by(name, lang, age, gender, direction) %>%
  summarise(fcr = mean(fcr, na.rm = T))

kids_large_table <- kids_spread %>%
  left_join(kids_fcr_spread)

kids_large_table %>%
  write_csv("large_table_kids.csv")

Visualizing Reversal Effect

I want to try to visualize reversal effects a different way. Maybe this.

# Get participant-level data
kids_fcr2 <- kids_fcr %>%
  group_by(name, age, lang, direction) %>%
  summarise(fcr = mean(fcr))

# reversal_effect_lm <- lmer(fcr ~ age + lang * direction + (1|name), data = kids_fcr2)
# summary(reversal_effect_lm)

ggplot(kids_fcr2, aes(x = direction, y = fcr, color = lang, fill = lang)) +
  geom_point() +
  geom_line(aes(group = name)) +
  facet_grid(. ~ lang) + 
  scale_y_continuous(limits = c(-1,1)) +
  theme_bw()

Or a reversal effect chart? Okay, so this chart tells us overall there really wasn’t much of a reversal effect for SE babies, they’re all hovering around 0. Interesting. While there seems to be a reversal effect for NSE babies where they look at the face more during reversed stories!

# Get participant-level data
kids_fcr3 <- kids_fcr2 %>%
  spread(direction, fcr) %>%
  group_by(name, age, lang) %>%
  mutate(diff = forward - reversed)

ggplot(kids_fcr3, aes(x = age, y = diff, color = lang)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  scale_y_continuous(limits = c(-1,1)) +
  theme_bw() +
  ggtitle("Reversal Effect") +
  ylab("Forward FCR - Reversed FCR")

And within-subjects variation here:

# First get the mean of each trial, THEN the participant-level means
within_subjects <- kids_fcr %>%
  group_by(name, lang, direction, story, repetition) %>%
  summarise(fcr = mean(fcr, na.rm = TRUE),
            count = n()) %>%
  group_by(name, lang, direction) %>%
  summarise(mean = mean(fcr, na.rm = TRUE),
            se = sd(fcr, na.rm = TRUE)/sqrt(n()),
            count = n())
# Then spread out mean and SE columns by direction
within_subjects_means <- within_subjects %>%
  select(-se, -count) %>%
  spread(direction, mean, sep = "_")
within_subjects_se <- within_subjects %>%
  select(-mean, -count) %>%
  spread(direction, se, sep = "SE")
within_subjects <- left_join(within_subjects_means, within_subjects_se, by = c("name","lang"))

# Now let's plot
lims <- c(-1,1)
within_subjects %>%
  ggplot(aes(x = direction_forward, y = direction_reversed, color = lang)) +
    geom_abline() +
  geom_point(size = 2) + 
  geom_errorbar(aes(ymin=direction_reversed-directionSEreversed, ymax=direction_reversed+directionSEreversed)) +
  geom_errorbarh(aes(xmin=direction_forward-directionSEforward, xmax=direction_forward+directionSEforward)) +
  theme_bw() +
  theme(aspect.ratio = 1) +
  scale_x_continuous("forward", limits = c(-1,1)) +
  scale_y_continuous("reversed", limits = c(-1,1)) +
  ggtitle("FCR Means") +
  facet_wrap("lang")

And a classic box/error plot with age collapsed.

kids_fcr2 %>%
  group_by(lang, direction) %>%
  summarise(fcr_mean = mean(fcr),
            sd = sd(fcr),
            n = n(),
            se = sd/sqrt(n)) %>%
  ggplot(aes(x = lang, y = fcr_mean, fill = direction)) +
  geom_bar(stat = "identity", position = position_dodge()) + 
  geom_errorbar(aes(ymin = fcr_mean-se, ymax = fcr_mean+se), position = position_dodge(0.9), width = 0.2) +
  scale_y_continuous(limits = c(-0.5, 0.5)) +
  theme_linedraw()

library(extrafont)

Registering fonts with R

# For making the babies/adults charts: 
kids_fcr2 %>%
    add_column(group = 'children') %>%
    write_csv("fcr_individual_values_children.csv")

kids_fcr2 %>%
  group_by(lang, direction) %>%
  summarise(fcr_mean = mean(fcr),
            sd = sd(fcr),
            n = n(),
            se = sd/sqrt(n)) %>%
  add_column(group = 'children') %>%
  write_csv("fcr_chart_children.csv")


kids_fcr2 %>%
  group_by(lang, direction) %>%
  summarise(fcr_mean = mean(fcr),
            sd = sd(fcr),
            n = n(),
            se = sd/sqrt(n)) %>%
  ggplot(aes(x = lang, y = fcr_mean, color = direction, fill = direction, group = direction)) +
  geom_hline(yintercept = 0, size = 0.5) +
  geom_point(size = 6, position = position_dodge(width = 0.4)) +
  geom_errorbar(aes(ymin = fcr_mean-se, ymax = fcr_mean+se), 
                size = 2, 
                position = position_dodge(0.4), 
                width = 0.3) +
  scale_y_continuous(limits = c(-0.5, 0.5)) +
  theme_linedraw() +
  theme(text = element_text(size = 30),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        axis.title.y = element_blank(),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        panel.border = element_rect(size = 2),
        axis.ticks.y = element_line(size = 0.5),
        panel.grid.major.y = element_line(size = 0.5, color = "light gray", linetype = "dashed")) +
  guides(color = FALSE, fill = FALSE)

Heat Maps

And now heat maps!

heatmap_kids <- kids %>%
  filter(aoi %in% midline) %>%
  ungroup() %>%
  group_by(lang, name, direction, aoi) %>%
  summarise(percent = mean(percent, na.rm=TRUE)) %>%
  group_by(lang, direction, aoi) %>%
  summarise(percent = mean(percent, na.rm=TRUE)) %>%
  ungroup() %>%
  mutate(aoi = factor(aoi, levels = c("Belly","BelowChest","MidChestBottom","MidChestCenter","MidChestTop",
            "MidFaceBottom","MidFaceCenter","MidFaceTop")))

ggplot(heatmap_kids, aes(x = lang, y = aoi)) +
  geom_tile(aes(fill=percent),color="lightgray",na.rm=TRUE) + 
#  scale_fill_viridis(option = "viridis", direction=-1, limits = c(0,.7), labels = percent, name = "looking time") +
    scale_fill_gradient(low = "#ffffff", high = "#08519c", space = "Lab", limits = c(0,.52), labels = percent, name = "looking time", na.value = "grey50") +
  theme_bw() +
  theme(strip.text.x = element_text(size = 11, color = "black", face = "italic"), 
        strip.background = element_rect(colour = "white", fill = "white"),
        panel.grid.major = element_line(color = "white")) +
  facet_grid(. ~ direction) +
  ylab("") + xlab("") + ggtitle("Eye Gaze Heat Map, by Direction") + 
  scale_y_discrete(expand=c(0,0)) +
  scale_x_discrete(expand = c(0,0))

ggplot(heatmap_kids, aes(x = direction, y = aoi)) +
  geom_tile(aes(fill=percent),color="lightgray",na.rm=TRUE) + 
#  scale_fill_viridis(option = "viridis", direction=-1, limits = c(0,.7), labels = percent, name = "looking time") +
    scale_fill_gradient(low = "#ffffff", high = "#08519c", space = "Lab", limits = c(0,.52), labels = percent, name = "looking time", na.value = "grey50") +
  theme_bw() +
  theme(strip.text.x = element_text(size = 11, color = "black", face = "italic"), 
        strip.background = element_rect(colour = "white", fill = "white"),
        panel.grid.major = element_line(color = "white")) +
  facet_grid(. ~ lang) +
  ylab("") + xlab("") + ggtitle("Eye Gaze Heat Map, by Language") + 
  scale_y_discrete(expand=c(0,0)) +
  scale_x_discrete(expand = c(0,0))

Collaped by direction (new)

heatmap_kids2 <- kids %>%
  filter(aoi %in% midline) %>%
  ungroup() %>%
  group_by(lang, name, aoi) %>%
  summarise(percent = mean(percent, na.rm=TRUE)) %>%
  group_by(lang, aoi) %>%
  summarise(percent = mean(percent, na.rm=TRUE)) %>%
  ungroup() %>%
  mutate(aoi = factor(aoi, levels = c("Belly","BelowChest","MidChestBottom","MidChestCenter","MidChestTop",
            "MidFaceBottom","MidFaceCenter","MidFaceTop")))

ggplot(heatmap_kids2, aes(x = lang, y = aoi)) +
  geom_tile(aes(fill=percent),color="lightgray",na.rm=TRUE) + 
#  scale_fill_viridis(option = "viridis", direction=-1, limits = c(0,.7), labels = percent, name = "looking time") +
    scale_fill_gradient(low = "#ffffff", high = "#08519c", space = "Lab", limits = c(0,.52), labels = percent, name = "looking time", na.value = "grey50") +
  theme_bw() +
  theme(strip.text.x = element_text(size = 11, color = "black", face = "italic"), 
        strip.background = element_rect(colour = "white", fill = "white"),
        panel.grid.major = element_line(color = "white")) +
  ylab("") + xlab("") + ggtitle("Eye Gaze Heat Map, by Language (Collapsed by Direction") + 
  scale_y_discrete(expand=c(0,0)) +
  scale_x_discrete(expand = c(0,0))

For poster?

ggplot(heatmap_kids, aes(x = direction, y = aoi)) +
  geom_tile(aes(fill=percent), 
            color="dark gray", 
            size = 0.25, 
            na.rm=T, 
            height = rep(c(10,4,1,1,1,1,1,1),4)
            ) + 
  scale_fill_gradient(low = "#ffffff", 
                      high = "#08519c", 
                      space = "Lab", 
                      limits = c(0,.52), 
                      labels = percent, 
                      name = "looking time", 
                      na.value = "grey50") +
  facet_grid(. ~ lang) +
  ylab("") + xlab("") + ggtitle("Eye Gaze Heat Map, by Language") + 
  scale_y_discrete(expand=c(0,0)) +
  scale_x_discrete(expand=c(0,0)) +
  theme_bw() +
  theme(text = element_text(size = 20),
        axis.title.y = element_blank(),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.x = element_blank(),
        axis.ticks.y = element_blank(),
        strip.text = element_blank(),
        panel.border = element_rect(size = 2),
        title = element_blank()) + 
  guides(color = FALSE, fill = FALSE)

# All Data
#Here's all AOI data. 

kids %>%
  ungroup() %>%
  group_by(lang, name, direction, aoi) %>%
  summarise(percent = mean(percent, na.rm=TRUE)) %>%
  group_by(lang, direction, aoi) %>%
  summarise(percent = mean(percent, na.rm=TRUE)) %>%
  openxlsx::write.xlsx("~/Desktop/kids_by_direction.xlsx")

Note: zip::zip() is deprecated, please use zip::zipr() instead

kids %>%
  ungroup() %>%
  group_by(lang, name, direction, aoi) %>%
  summarise(percent = mean(percent, na.rm=TRUE)) %>%
  group_by(lang, direction, aoi) %>%
  summarise(percent = mean(percent, na.rm=TRUE)) %>%
  group_by(lang, aoi) %>%
  summarise(percent = mean(percent, na.rm=TRUE)) %>%
  openxlsx::write.xlsx("~/Desktop/kids_collapsed.xlsx")

Discussion

No big changes from the ICSLA abstract. Good!

The interpretation here is that:

• All kids looked equally at all videos regardless of language or direction. Age did have an effect so we used age in our models. Good!
• SE kids continue to be strong face-lookers compared to NSE kids. (Same as ICSLA)
• There is no reversal effect. (Same as ICSLA)

That doesn’t mean both groups of children don’t care about reversal. On the contrary. We can hypothesize that SE kids have efficient gaze behavior and are resilient to reversal; while NSE kids already are “inefficient” and changing the video stimulus isn’t going to help. But how do we test that? Maybe let’s look at within-subject variation.

Correlations

Let’s try correlations.

# Let's try correlations
kids_nse <- kids %>% 
  filter(aoi %in% midline) %>%
  filter(lang == "NSE") %>%
  group_by(name, direction, aoi) %>% 
  summarise(percent = mean(percent)) %>%
  ungroup() %>%
  mutate(direction = case_when(
    direction == "forward" ~ "fw",
    direction == "reversed" ~ "rv"
  )) %>% 
  unite(aoi2, direction, aoi, sep = "_") %>%
  spread(aoi2, percent) %>%
  select(-name)

kids_se <- kids %>% 
  filter(aoi %in% midline) %>%
  filter(lang == "SE") %>%
  group_by(name, direction, aoi) %>% 
  summarise(percent = mean(percent)) %>%
  ungroup() %>%
  mutate(direction = case_when(
    direction == "forward" ~ "fw",
    direction == "reversed" ~ "rv"
  )) %>% 
  unite(aoi2, direction, aoi, sep = "_") %>%
  spread(aoi2, percent) %>%
  select(-name)

ggcorr(kids_nse, label = TRUE, label_size = 5, label_round = 2, label_alpha = TRUE, hjust = 0.9, size = 5, color = "grey50", layout.exp = 1) + ggtitle("NSE")

ggcorr(kids_se, label = TRUE, label_size = 5, label_round = 2, label_alpha = TRUE, hjust = 0.9, size = 5, color = "grey50", layout.exp = 1) + ggtitle("SE")

library(corrr)
kids_nse %>% correlate() %>% network_plot(min_cor=0.6) + ggtitle("NSE Children")

Correlation method: 'pearson'
Missing treated using: 'pairwise.complete.obs'

kids_se %>% correlate() %>% network_plot(min_cor=0.6) + ggtitle("SE Children")

Correlation method: 'pearson'
Missing treated using: 'pairwise.complete.obs'

XY Space Data

We’ll load the data from the childxydata.feather file made in 06rawxydata.Rmd. So any new kids, please run the first code block in 06 to include it. Then we’ll keep all the kids we also have in the AOI data group.

included <- kids %>%
  ungroup() %>%
  select(name) %>% 
  distinct() %>%
  unlist()

xydata <- read_feather("../Child Data/childxydata.feather") %>%
  rename(name = participant) %>%
  filter(name %in% included)

# Get ages
ages <- read_csv("childrenages.csv") %>%
  rename(name = participant)
xydata <- xydata %>% left_join(ages, by = "name") %>%
  mutate(age = age*12) %>%
  mutate(agegroup = case_when(
    age <= 8.99 ~ "younger",
    age >= 9.0 & age < 15 ~ "older"
  )) %>%
  mutate(language = case_when(
    language == "EnglishExposed" ~ "NSE",
    language == "SignLanguageExposed" ~ "SE"
  )) %>%
  rename(lang = language) %>%
  select(name, group, gender, lang, condition, mark, trial, repetition, x, y, age, agegroup) %>%
  separate(condition, into = c("story", "clip", "direction")) %>%
  unite("story", c("story", "clip")) %>%
  mutate(direction = case_when(
    direction == "ER" ~ "reversed",
    direction == "FW" ~ "forward"
  )) %>%
  mutate(name = factor(name),
         group = factor(group),
         gender = factor(gender),
         lang = factor(lang),
         story = factor(story),
         direction = factor(direction),
         mark = factor(mark),
         trial = factor(trial),
         repetition = factor(repetition),
         agegroup = factor(agegroup))

Overall Looking

Let’s check that we have no significant group or condition differences in terms of valid (not empty) data points collected. This is same as “Global Looking” we have above, really, but w raw xy data.

xy_overall <- xydata %>%
  filter(!is.na(x)) %>%
  group_by(name, age, lang, direction, story, repetition) %>%
  summarise(data_points = n()) # gets total looking percent for each trial for each baby

# Table of means
xy_overall %>% 
  group_by(name, lang, direction) %>%
  summarise(data_points = mean(data_points)) %>% # get average looking percent for each baby
  group_by(lang, direction) %>%
  summarise(mean_data_points = mean(data_points),
            count = n(),
            sd = sd(data_points),
            se = sd/sqrt(count)) %>%
  select(-sd) %>%
  print()

ggplot(xy_overall, aes(x = age, y = data_points, color = direction, fill = direction)) + 
  geom_jitter(alpha = 0.5) +
  facet_grid(. ~ lang) +
  geom_smooth(method = "lm", se = FALSE) +
  ggtitle("Data Points") +
  xlab("age (months)") +
  ylab("data points recorded") + 
  theme_bw() 

Description.

overall_xy_lm <- lmer(data_points ~ age + lang * direction + (direction|name) + (direction|story), data = xy_overall)
summary(overall_xy_lm) 

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: data_points ~ age + lang * direction + (direction | name) + (direction |  
    story)
   Data: xy_overall

REML criterion at convergence: 7516.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.0133 -0.5789  0.1903  0.7322  2.2752 

Random effects:
 Groups   Name              Variance  Std.Dev. Corr 
 name     (Intercept)       1.888e+04 137.3995      
          directionreversed 1.995e-01   0.4466 0.97 
 story    (Intercept)       4.162e+04 204.0169      
          directionreversed 9.655e+03  98.2577 -0.44
 Residual                   6.768e+04 260.1538      
Number of obs: 535, groups:  name, 35; story, 8

Fixed effects:
                         Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)               524.679    117.283  28.993   4.474 0.000109 ***
age                         2.476      1.382  32.811   1.791 0.082539 .  
langSE                      4.408     57.810  32.196   0.076 0.939689    
directionreversed         -17.953     49.870  12.327  -0.360 0.724939    
langSE:directionreversed   10.149     47.247 487.388   0.215 0.830003    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) age    langSE drctnr
age         -0.694                     
langSE      -0.269 -0.024              
dirctnrvrsd -0.296 -0.002  0.227       
lngSE:drctn  0.114  0.005 -0.405 -0.553

confint(overall_xy_lm)

Computing profile confidence intervals ...
non-monotonic profile for .sig05bad spline fit for .sig05: falling back to linear interpolation

                                2.5 %      97.5 %
.sig01                     97.5119565 179.9039190
.sig02                     -1.0000000   1.0000000
.sig03                      0.0000000  48.5054545
.sig04                    121.2560632 335.9687619
.sig05                     -0.8443818   0.5708122
.sig06                     28.1071476 180.6104114
.sigma                    244.4048523 277.1006287
(Intercept)               305.0235868 745.6515276
age                        -0.1729330   5.0838822
langSE                   -105.1902957 115.6896029
directionreversed        -110.6463260  75.2664550
langSE:directionreversed  -82.8656734 102.4687275

#ggcoef(overall_xy_lm)

XY Data LMMs

Now we’re going to run LMMs on babies’ raw:

• horizontal spread (middle 50% of x data; xIQR)
• vertical spread (middle 50% of y data; yIQR)
• viewing area (A = middle-x * middle-y; area)

But to do this we first trim each kid’s data, getting rid of the first 60 samples (0.50 secs) of each trial.

xydata <- xydata %>%
  group_by(name,trial) %>%
  slice(30:n())

iqr <- xydata %>%
  group_by(name, age, lang, story, direction, trial) %>%
  summarise(xIQR = IQR(x,na.rm=TRUE),
                   yIQR = IQR(y,na.rm=TRUE),
                   xmed = median(x, na.rm=TRUE),
                   ymed = median(y, na.rm=TRUE),
                   area = xIQR*yIQR)
head(iqr,20)

NA

Middle X

Description.

xiqr_mean <- iqr %>% 
  group_by(lang, direction, name) %>%
  summarise(xIQR = mean(xIQR, na.rm = T)) %>% # gets average looking percent for each baby
  group_by(lang, direction) %>%
  summarise(mean_xIQR = mean(xIQR), # gets group averages
            count = n(),
            sd = sd(xIQR),
            se = sd/sqrt(count)) %>%
  select(-sd) %>%
  print()

# Plot
ggplot(iqr, aes(x = age, y = xIQR, color = direction, fill = direction)) + 
  geom_jitter(alpha = 0.5) + 
  geom_smooth(method = "lm", se = FALSE) +
  theme_bw() + 
  theme(panel.grid.major.x = element_blank()) +
  facet_grid(. ~ lang) +
  ggtitle("Horizontal Spread") +
  xlab("") +
  ylab("xIQR")

ggplot(xiqr_mean, aes(x = lang, y = mean_xIQR, fill = direction)) +
  geom_bar(stat = "identity", position = position_dodge()) + 
  geom_errorbar(aes(ymin = mean_xIQR-se, ymax = mean_xIQR+se), position = position_dodge(0.9), width = 0.2) +
  theme_linedraw()

xiqr_lm <- lmer(xIQR ~ age + lang * direction + (1|name) + (1|story), data = iqr)
summary(xiqr_lm)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: xIQR ~ age + lang * direction + (1 | name) + (1 | story)
   Data: iqr

REML criterion at convergence: 5145.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.3266 -0.4336 -0.1667  0.1542 12.5365 

Random effects:
 Groups   Name        Variance Std.Dev.
 name     (Intercept)  61.846   7.864  
 story    (Intercept)   8.719   2.953  
 Residual             878.268  29.636  
Number of obs: 534, groups:  name, 35; story, 8

Fixed effects:
                          Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)               38.13628    6.85205  39.38839   5.566 2.02e-06 ***
age                       -0.06613    0.09760  32.39698  -0.678    0.503    
langSE                    -0.89305    4.55252  69.37559  -0.196    0.845    
directionreversed          3.17906    3.95964 434.80176   0.803    0.422    
langSE:directionreversed  -2.67056    5.25883 457.77145  -0.508    0.612    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) age    langSE drctnr
age         -0.851                     
langSE      -0.378 -0.003              
dirctnrvrsd -0.289 -0.002  0.440       
lngSE:drctn  0.214  0.007 -0.573 -0.759

confint(xiqr_lm)

Computing profile confidence intervals ...

                               2.5 %     97.5 %
.sig01                     3.5754861 11.2093528
.sig02                     0.0000000  7.2383649
.sigma                    27.8313919 31.5418940
(Intercept)               24.9210180 51.3220302
age                       -0.2546145  0.1226814
langSE                    -9.6578128  7.9002048
directionreversed         -4.5353782 10.9675126
langSE:directionreversed -13.0257923  7.5654776

#ggcoef(xiqr_lm)

Middle Y

Description.

yiqr_mean <- iqr %>% 
  group_by(lang, direction, name) %>%
  summarise(yIQR = mean(yIQR, na.rm = T)) %>% # gets average looking percent for each baby
  group_by(lang, direction) %>%
  summarise(mean_yIQR = mean(yIQR), # gets group averages
            count = n(),
            sd = sd(yIQR),
            se = sd/sqrt(count)) %>%
  select(-sd) %>%
  print()

# Plot
ggplot(iqr, aes(x = age, y = yIQR, color = direction, fill = direction)) + 
  geom_jitter(alpha = 0.5) + 
  geom_smooth(method = "lm", se = FALSE) +
  theme_bw() + 
  theme(panel.grid.major.x = element_blank()) +
  facet_grid(. ~ lang) +
  ggtitle("Vertical Spread") +
  xlab("") +
  ylab("yIQR")

ggplot(yiqr_mean, aes(x = lang, y = mean_yIQR, fill = direction)) +
  geom_bar(stat = "identity", position = position_dodge()) + 
  geom_errorbar(aes(ymin = mean_yIQR-se, ymax = mean_yIQR+se), position = position_dodge(0.9), width = 0.2) +
  theme_linedraw()

yiqr_lm <- lmer(yIQR ~ age + lang * direction + (1|name) + (1|story), data = iqr)
summary(yiqr_lm)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: yIQR ~ age + lang * direction + (1 | name) + (1 | story)
   Data: iqr

REML criterion at convergence: 5653.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.5584 -0.5076 -0.2419  0.1108  5.9858 

Random effects:
 Groups   Name        Variance Std.Dev.
 name     (Intercept)  300.035 17.322  
 story    (Intercept)    9.853  3.139  
 Residual             2251.577 47.451  
Number of obs: 534, groups:  name, 35; story, 8

Fixed effects:
                         Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)               59.4257    12.9691  37.0247   4.582 5.08e-05 ***
age                       -0.1215     0.1887  32.7752  -0.644   0.5242    
langSE                   -10.4519     8.3316  56.6168  -1.254   0.2148    
directionreversed         10.9450     6.2897 400.4069   1.740   0.0826 .  
langSE:directionreversed -10.0349     8.3652 432.1321  -1.200   0.2310    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) age    langSE drctnr
age         -0.870                     
langSE      -0.367  0.000              
dirctnrvrsd -0.243 -0.001  0.381       
lngSE:drctn  0.179  0.006 -0.498 -0.755

confint(yiqr_lm)

Computing profile confidence intervals ...

                               2.5 %     97.5 %
.sig01                    10.8614437 23.1060826
.sig02                     0.0000000  9.9939596
.sigma                    44.5593590 50.4957715
(Intercept)               34.4378120 84.3887649
age                       -0.4854389  0.2425792
langSE                   -26.4843811  5.5982767
directionreversed         -1.3514878 23.2682377
langSE:directionreversed -26.4541999  6.2990337

#ggcoef(yiqr_lm)

Viewing Area

Description.

area_mean <- iqr %>% 
  group_by(lang, direction, name) %>%
  summarise(area = mean(area, na.rm = T)) %>% # gets average looking percent for each baby
  group_by(lang, direction) %>%
  summarise(area_mean = mean(area), # gets group averages
            count = n(),
            sd = sd(area),
            se = sd/sqrt(count)) %>%
  select(-sd) %>%
  print()

# Plot
ggplot(iqr, aes(x = age, y = area, color = direction, fill = direction)) + 
  geom_jitter(alpha = 0.5) + 
  geom_smooth(method = "lm", se = FALSE) +
  theme_bw() + 
  theme(panel.grid.major.x = element_blank()) +
  facet_grid(. ~ lang) +
  ggtitle("Viewing Area") +
  xlab("") +
  ylab("Area (px^2)")

ggplot(area_mean, aes(x = lang, y = area_mean, fill = direction)) +
  geom_bar(stat = "identity", position = position_dodge()) + 
  geom_errorbar(aes(ymin = area_mean-se, ymax = area_mean+se), position = position_dodge(0.9), width = 0.2) +
  theme_linedraw()

area_lm <- lmer(area ~ age + lang * direction + (1|name) + (1|story), data = iqr)
summary(area_lm)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: area ~ age + lang * direction + (1 | name) + (1 | story)
   Data: iqr

REML criterion at convergence: 10618.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.0099 -0.2801 -0.1786 -0.0233 15.8524 

Random effects:
 Groups   Name        Variance Std.Dev.
 name     (Intercept)  1702595 1304.8  
 story    (Intercept)    31912  178.6  
 Residual             27599640 5253.5  
Number of obs: 534, groups:  name, 35; story, 8

Fixed effects:
                          Estimate Std. Error        df t value Pr(>|t|)  
(Intercept)               2605.197   1165.076    39.619   2.236    0.031 *
age                         -9.887     16.741    33.451  -0.591    0.559  
langSE                     181.489    785.881    73.046   0.231    0.818  
directionreversed          881.525    691.159   362.314   1.275    0.203  
langSE:directionreversed -1285.676    920.304   402.572  -1.397    0.163  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) age    langSE drctnr
age         -0.859                     
langSE      -0.383 -0.003              
dirctnrvrsd -0.297 -0.001  0.442       
lngSE:drctn  0.218  0.007 -0.581 -0.752

confint(area_lm)

Computing profile confidence intervals ...
Last two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepunexpected decrease in profile: using minstepunexpected decrease in profile: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepLast two rows have identical or NA .zeta values: using minstepnon-monotonic profile for .sig02bad spline fit for .sig02: falling back to linear interpolationcollapsing to unique 'x' values

                               2.5 %     97.5 %
.sig01                     522.19135 1876.55263
.sig02                       0.00000        Inf
.sigma                    4936.29250 5587.86945
(Intercept)                368.63639 4849.05524
age                        -42.13398   22.39641
langSE                   -1336.66470 1688.46326
directionreversed         -478.08996 2222.87739
langSE:directionreversed -3075.97943  522.77446

#ggcoef(area_lm)

Plotting Viewing Area

medians <- iqr %>%
  group_by(name,lang,direction) %>%
  summarise(xIQR = mean(xIQR,na.rm=TRUE),
                   yIQR = mean(yIQR,na.rm=TRUE),
                   xmed = mean(xmed,na.rm=TRUE),
                   ymed = mean(ymed,na.rm=TRUE)) %>%
  group_by(lang,direction) %>% 
  summarise(xIQR = mean(xIQR,na.rm=TRUE),
                   yIQR = mean(yIQR,na.rm=TRUE),
                   x = mean(xmed,na.rm=TRUE),
                   y = mean(ymed,na.rm=TRUE)) %>%
  mutate(y = y*-1,
         xmin = x-(xIQR/2),
         xmax = x+(xIQR/2),
         ymin = y-(yIQR/2),
         ymax = y+(yIQR/2))
img <- png::readPNG("cindy.png")
g <- grid::rasterGrob(img, interpolate=TRUE, width=unit(1,"npc"), height=unit(1,"npc")) 
ggplot(medians, aes(fill=direction,color=direction)) +
  annotation_custom(g, xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=Inf) +
  geom_rect(aes(xmin=xmin,ymin=ymin,xmax=xmax,ymax=ymax),alpha=.1) + 
  theme_linedraw() +
  scale_x_continuous(limits = c(0,1080), expand = c(0, 0)) +
  scale_y_continuous(limits = c(-720,0), expand = c(0, 0)) +
  facet_wrap("lang")

XY Space Data - Multiple Plots

First let’s prep the data.

multiples <- xydata %>%
  filter(!is.na(x)) %>%
  filter(!is.na(y)) %>%
  group_by(name, age, lang, story, direction, trial) %>%
  summarise(xIQR = IQR(x,na.rm=TRUE),
            yIQR = IQR(y,na.rm=TRUE),
            xmed = median(x, na.rm=TRUE),
            ymed = median(y, na.rm=TRUE),
            area = xIQR*yIQR,
            x_90 = quantile(x, .95, na.rm=TRUE) - quantile(x, .05, na.rm=TRUE),
            y_90 = quantile(y, .95, na.rm=TRUE) - quantile(y, .05, na.rm=TRUE),
            area_90 = (x_90) * (y_90),
            x_mean = mean(x, na.rm = TRUE),
            y_mean = mean(y, na.rm = TRUE),
            x_sd = sd(x, na.rm = TRUE),
            y_sd = sd(y, na.rm = TRUE),
            x_1sd = (x_mean+x_sd) - (x_mean-x_sd),
            y_1sd = (y_mean+y_sd) - (y_mean-y_sd),
            area_1sd = x_1sd * y_1sd,
            x_2sd = (x_mean+(x_sd*2)) - (x_mean-(x_sd*2)),
            y_2sd = (y_mean+(y_sd*2)) - (y_mean-(y_sd*2)),
            area_2sd = x_2sd * y_2sd) %>%
  group_by(name, lang, direction) %>%
  summarise_if(is.double, funs(mean), na.rm = T) %>%
  group_by(lang, direction) %>%
  summarise_if(is.double, funs(mean), na.rm = T)

funs() is soft deprecated as of dplyr 0.8.0
Please use a list of either functions or lambdas: 

  # Simple named list: 
  list(mean = mean, median = median)

  # Auto named with `tibble::lst()`: 
  tibble::lst(mean, median)

  # Using lambdas
  list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
This warning is displayed once per session.

img <- png::readPNG("cindy.png")
g <- grid::rasterGrob(img, interpolate=TRUE, width=unit(1,"npc"), height=unit(1,"npc")) 

IQR (Middle 50%)

Let’s see.

curr_data <- multiples %>% 
  ungroup() %>%
  select(lang, direction, xmed, ymed, xIQR, yIQR) %>%
  group_by(lang, direction) %>%
  summarise(xmin = xmed-(xIQR/2),
         xmax = xmed+(xIQR/2),
         ymin = -1*(ymed-(yIQR/2)),
         ymax = -1*(ymed+(yIQR/2)))

ggplot(curr_data, aes(fill=direction,color=direction)) +
  annotation_custom(g, xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=Inf) +
  geom_rect(aes(xmin=xmin,ymin=ymin,xmax=xmax,ymax=ymax),alpha=.1) + 
  theme_linedraw() +
  scale_x_continuous(limits = c(0,1080), expand = c(0, 0)) +
  scale_y_continuous(limits = c(-720,0), expand = c(0, 0)) +
  facet_wrap("lang")

Middle 90%

So I calculated the average median across, and the middle 90% of the data.

curr_data <- multiples %>% 
  ungroup() %>%
  select(lang, direction, xmed, ymed, x_90, y_90) %>%
  group_by(lang, direction) %>%
  summarise(xmin = xmed-(x_90/2),
         xmax = xmed+(x_90/2),
         ymin = -1*(ymed-(y_90/2)),
         ymax = -1*(ymed+(y_90/2)))

ggplot(curr_data, aes(fill=direction,color=direction)) +
  annotation_custom(g, xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=Inf) +
  geom_rect(aes(xmin=xmin,ymin=ymin,xmax=xmax,ymax=ymax),alpha=.1) + 
  theme_linedraw() +
  scale_x_continuous(limits = c(0,1080), expand = c(0, 0)) +
  scale_y_continuous(limits = c(-720,0), expand = c(0, 0)) +
  facet_wrap("lang")

# ggplot(filter(curr_data, lang == "NSE"), aes(fill=direction,color=direction)) +
#   annotation_custom(g, xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=Inf) +
#   geom_rect(aes(xmin=xmin,ymin=ymin,xmax=xmax,ymax=ymax),alpha=.2, size = 1) +
#   theme_linedraw() +
#   scale_x_continuous(limits = c(0,1080), expand = c(0, 0)) +
#   scale_y_continuous(limits = c(-720,0), expand = c(0, 0))
# 
# 
# ggplot(filter(curr_data, lang == "SE"), aes(fill=direction,color=direction)) +
#   annotation_custom(g, xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=Inf) +
#   geom_rect(aes(xmin=xmin,ymin=ymin,xmax=xmax,ymax=ymax),alpha=.2, size = 1) +
#   theme_linedraw() +
#   scale_x_continuous(limits = c(0,1080), expand = c(0, 0)) +
#   scale_y_continuous(limits = c(-720,0), expand = c(0, 0))

±1 SD (Middle 68%)

So this is using the mean of the means, plus or minus one SD. This is equivalent to middle 68%.

curr_data <- multiples %>% 
  ungroup() %>%
  select(lang, direction, x_mean, y_mean, x_1sd, y_1sd) %>%
  group_by(lang, direction) %>%
  summarise(xmin = x_mean-(x_1sd/2),
         xmax = x_mean+(x_1sd/2),
         ymin = -1*(y_mean-(y_1sd/2)),
         ymax = -1*(y_mean+(y_1sd/2)))

ggplot(curr_data, aes(fill=direction,color=direction)) +
  annotation_custom(g, xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=Inf) +
  geom_rect(aes(xmin=xmin,ymin=ymin,xmax=xmax,ymax=ymax),alpha=.1) + 
  theme_linedraw() +
  scale_x_continuous(limits = c(0,1080), expand = c(0, 0)) +
  scale_y_continuous(limits = c(-720,0), expand = c(0, 0)) +
  facet_wrap("lang")

±2 SD (Middle 96%)

And this is using the mean of the means, plus or minus two SD. This is equivalent to middle 96%.

curr_data <- multiples %>% 
  ungroup() %>%
  select(lang, direction, x_mean, y_mean, x_2sd, y_2sd) %>%
  group_by(lang, direction) %>%
  summarise(xmin = x_mean-(x_2sd/2),
         xmax = x_mean+(x_2sd/2),
         ymin = -1*(y_mean-(y_2sd/2)),
         ymax = -1*(y_mean+(y_2sd/2)))

ggplot(curr_data, aes(fill=direction,color=direction)) +
  annotation_custom(g, xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=Inf) +
  geom_rect(aes(xmin=xmin,ymin=ymin,xmax=xmax,ymax=ymax),alpha=.1) + 
  theme_linedraw() +
  scale_x_continuous(limits = c(0,1080), expand = c(0, 0)) +
  scale_y_continuous(limits = c(-720,0), expand = c(0, 0)) +
  facet_wrap("lang")